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We develop an all-electron path integral Monte Carlo (PIMC) method with free-particle nodes for 
warm dense matter and apply it to water and carbon plasmas. We thereby extend PIMC studies 
beyond hydrogen and helium to elements with core electrons. PIMC pressures, internal energies, and 
pair-correlation functions compare well with density functional theory molecular dynamics (DFT- 
MD) at temperatures of (2.5-7.5) x 10^ K and both methods together form a coherent equation of 
state (EOS) over a density-temperature range of 3-12 g/cm'' and 10^-10^ K. 
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The development of first-principles methodology for 
warm, dense matter (WDM) is one of the great challenges 
of modern materials theory. A need for rigorous sim- 
ulation of WDM has escalated with intensified interest 
in advanced energy technologies [li, physics and chem- 
istry of solar and extrasolar planets [4I , shock compressed 
matterlsl , and different types of plasma-matter interac- 
tions y. The standard first-principles method, Kohn- 
Sham density functional theory molecular dynamics [5J 
(DFT-MD), produces accurate equations of state in the 
lower temperature range of the WDM regime. The maxi- 
mum accessible temperature is limited, however, because 
the number of partially occupied orbitals eventually be- 
comes computationally intractable [6|. On the other 
hand, the many-body path integral Monte Carlo (PIMC) 
method [7| is naturally formulated to study high temper- 
ature dependence of materials. Ideally, PIMC and DFT 
together can produce a coherent equation of state for the 
entire WDM regime and cross- validate each other at com- 
monly accessible temperatures. However, PIMC has not 
yet been developed to study systems with core electrons. 
Indeed, PIMC studies up to now have been limited to 
plasma states of hydrogen |8l4lC| and helium ll|, |l2| . In 
this letter, we develop an all-electron PIMC method for 
first-row elements and combine results with DFT-MD to 
produce comprehensive equations of state for water and 
carbon in the WDM regime for a density-temperature 
range of 3-12 g/cm^ and 10^-10^ K. 

The central characteristic of a material in the WDM 
regime is that the electron-ion interaction becomes com- 
parable to the electron kinetic energy and, therefore, ef- 
fects of bonding, ionization, exchange-correlation, and 
quantum degeneracy are all important. The analytic 
methods of condensed matter and plasma physics [13| are 
typically not reliable without experimental input. One 
must turn to the the numerical, first-principles PIMC and 
DFT-MD methods which accurately capture the many- 
body physics in the WDM regime without empirical pa- 
rameters or corrections. However, first-principles meth- 
ods utilize certain approximations and one must compare 
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FIG. 1: (color online) Comparison of excess pressure relative 
to the ideal Fermi gas plotted as a function of temperature 
for water. 



with experimental data if available. 

The key approximation in DFT is that of the exchange- 
correlation potential, which describes all the many body 
interactions. The exchange-correlation potentials used in 
nearly all condensed matter calculations are constructed 
from zero temperature quantum Monte Carlo calcula- 
tions of the electron gas [1J| . In the WDM regime, tem- 
peratures are at or above the Fermi temperature and elec- 
trons are excited relative to their ground state. There- 
fore, without further validation, the exchange-correlation 
potential cannot be assumed to provide an accurate de- 
scription in the WDM regime. 

In DFT calculations it is also common to replace the 
core electrons in each atom with a pseudopotential. Typi- 
cally, highest accuracy is obtained with a non-local pseu- 
dopotential which depends on the energy and angular 
momentum components in core states. However, in the 
WDM regime, it is possible to excite electrons out of core 
levels. The pseudopotential approximation may break 



down and should always be compared with all-clcctron 
calculations. Additionally, finite-temperature DFT uses 
a Fermi-Dirac function to allow for thermal occupation 
of single-particle electronic states [l^, but requires an 
increasing number of bands with temperature, crippling 
its efficiency at extreme temperatures. Orbital-free den- 
sity functional methods aim to overcome such thermal 
band limitations, but several challenges remain to be 
solved (H. 

The PIMC method avoids the band structure calcula- 
tion and exchange-correlation approximation by being di- 
rectly defined from the path integral formulation of quan- 
tum statistics. PIMC stochastically solves the full finite- 
temperature quantum many-body problem by treating 
electrons and nuclei equally as paths and addresses all 
of WDM characteristics on an equal footing. All finite- 
temperature properties of a material are then readily cal- 
culated from the thermal density matrix. In contrast to 
DFT, PIMC efficiency increases with increasing temper- 
ature as particles become more classical and fewer time 
slices are needed, scaling inversely with temperature. A 
non-local pseudopotential formulation of PIMC docs not 
yet exist [l7| and this is why PIMC calculations so far 
have been limited to hydrogen and helium. PIMC calcu- 
lations presented here treat all electrons explicitly. 

The only uncontrolled approximation in PIMC is that 
of the nodal surface to deal with the fermion sign prob- 
lem. Unchecked, the fermion sign problem leads to a 
cancellation of positive and negative contributions to the 
density matrix which causes large fluctuations in com- 
puted averages. One solution to this problem is the so- 
called fixed-node approximation in which the location of 
the nodes are fixed to a known trial nodal structure in 
order to guarantee positive contributions to the thermal 
density matrix. The form of the density matrix does not 
enter the PIMC computation, only the location of the 
nodes. 

The PIMC method we present here employs a free- 
particle nodal structure, which is expected to be accurate 
for systems that are almost fully ionized. One could as- 
sume accurate calculations of heavier elements requires 
very high temperatures where atomic cores are ionized 
also. However, for hydrogen, PIMC with free-particle 
nodes has provided reliable results at much lower tem- 
peratures where only partial ionization occurs [lO[. The 
PIMC results presented for water and carbon here will 
demonstrate that accurate pressures and energies are ob- 
tained for temperatures so low that the Is states are fully 
occupied and the 2s states are partially occupied. Anal- 
yses of the DFT band occupations show that as the 2s 
states become less than 60% occupied for T > 2.5 x 10^ 
K, PIMC and DFT results agree. 

In order to explain this result, we first note that no 
nodes are needed to describe an isolated, doubly occu- 
pied Is state. Our results for water and carbon indicate 
that free particle nodes also work in cases where the Is 
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FIG. 2: (color online) Comparison of excess pressures relative 
to the ideal Fermi gas plotted as a function of temperature 
for carbon. 



state is doubly occupied and all other electrons are ion- 
ized. This may be because only one orbital out of many in 
the Slater determinant is not characterized well. As the 
occupation of the 2s state exceeds 60% at lower temper- 
atures, the PIMC pressures and energies become inaccu- 
rate because free particle nodes cannot yield the correct 
shell structure around the nucleus [l8|. Our results will 
show, for first-row elements, free particle nodes remain 
sufficiently accurate at low enough temperatures to over- 
lap with the highest temperature DFT data. This allows 
the two methods to cross- validate each other and form a 
single coherent equation of state for all temperatures. 

As a first application to test our method, we study 
water because it is one of the most prevalent materials 
in nature and knowledge of its electronic properties in 
the WDM regime is crucial for understanding aspects 
of astrophysical objects, such as the interiors of giant 
gas planets. Reports suggest Uranus, Neptune, Jupiter, 
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and Saturn contain significant amounts of water 
In addition to its rich solid and fiuid phases, water is 
known for its superionic and plasma phases as well as 
an insulator-to-metal transition at extreme densities and 
temperatures. Recent DFT-MD simulations [22| have 
computed the equation of state of water up to 2x10^ 
K and 15 g/cm'^, improving upon the older SESAME 
7150 [23[ table comprised of a number of analytic models 
and MD using empirical potentials. 

As a second application, we study carbon at high 
pressures and temperatures for its importance in future 
energy technologies. In inertial confined fusion experi- 
ments, carbon is used as an ablator for target capsules. 
The performance of the ablator is heavily dependent on 
the equation of state in the WDM regime [2J|. There 
have been a number of attempts to construct carbon 
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FIG. 3: (color online) Comparison of excess internal energies 
relative to the ideal Fermi gas plotted as a function of tem- 
perature for (a) carbon and (b) water. The lower density data 
has been shifted by a constant, -2, in both cases. 



equations of state in the WDM regime, including free 
energy models ^^ and DFT-MD ^, but they ulti- 
mately resort to more approximate Thomas-Fermi-based 
models that cannot describe any chemical bond. 

For our PIMC simulations, the Coulomb interaction is 
incorporated via pair density matrices derived from the 
eigenstates of the two-body Coulomb problem. A suffi- 
ciently small time step is determined by converging total 
energy as a function of time step until the energy changes 
less than 0.2%. For both water and carbon, we use a time 
step of 0.0078125 Ha^^ for temperatures below 5x10^ K 
and, for higher temperatures, the time step decreases as 
l/T while keeping at least four time slices in the path 
integral. In order to minimize finite size errors, the total 
energy is converged to better than 0.2% for a 24 atom 
simple cubic cell. 

The DFT-MD simulations were performed with ei- 
ther the ABINIT code [28| using all-electron pseudopo- 
tentials or with the Vienna ab initio simulation pack- 
age (VASP) [29[ using the projector augmented- wave 
method [30[. MD uses a NVT ensemble regulated with 
a Nose-Hoover thermostat. Exchange-correlation effects 
are described using the Perdew-Burke-Ernzerhof general- 
ized gradient approximation |3l| . Electronic wave func- 
tions are expanded in a plane-wave basis with a energy 
cutoff of at least 1500 eV for water and at least 900 eV 
for carbon in order to converge total energy to chemical 
accuracy. Size convergence tests indicate that total en- 
ergies are converged to better than 0.2% in a 24 atom 
simple cubic cell. Convergence of the number of orbitals 
to a prescribed thermal occupation of less than 1x10^^ 
requires up to 1500 bands at 7.5x10^ K for a 24-atom 
cell. All simulations are performed at the gamma-point 



FIG. 4: (color online) Nuclear pair-correlation functions for 
(a) carbon and (b) water. 



of the Brillouin zone, which converges total energy to 
better than 0.1% at relevant high temperatures. 

Figures 1 and 2 compare pressures obtained for wa- 
ter and carbon, respectively, from PIMC, DFT-MD, and 
Debye-Hiickel [32| simulations. Water is studied at fixed 
densities of 3.18 and 11.18 g/cm'^ and carbon is stud- 
ied at 4.17 g/cm'^ and 12.64 g/cm'^. The two densities 
in each case correspond to a pressure of 1 Mbar and 50 
Mbar at zero temperature. Pressures, P, are plotted rel- 
ative to a fully ionized Fermi gas of electrons and ions 
with pressure, Pq, in order to compare only the pres- 
sure contributions that result only from particle interac- 
tions. PIMC and DFT-MD rcsuhs for {P-Po)/Po agree 
to better than 0.03 in the range of 2.5x10^ to 7.5x10^ 
K. Convergence tests show that results are equally well 
converge in 24-atoin and 8-atom simulation cells. The 
excellent agreement allows for cross-validation which im- 
plies the zero temperature DFT exchange-correlation po- 
tential remains valid at high temperatures and that the 
free-particle nodal approximation is valid in PIMC when 
atoms are only partially ionized. The two methods have 
comparable computational cost in the overlap region, but 
DFT computational cost starts to become prohibitive be- 
yond 7.5x10^ K, and free particle nodes break down be- 
low 2.5x10^ K. 

In addition, Fig. 2 compares the instantaneous pres- 
sures obtained for a fixed configuration of carbon nu- 
clei at various electronic temperatures using PIMC, DFT 
with all electron pseudopotentials, and DFT with VASP 
PAW pseudopotentials. Agreement between PIMC and 
DFT with all electron pseudopotentials is very good from 
1x10^ to 2x10^ K. However, beyond 7.5x10^ K, PAW 
DFT no longer predicts the correct temperature depen- 
dence, indicating that the missing contributions of core 
excitations to the total energy become significant. All 



electron DFT is too computationally demanding to per- 
form calculations with moving nuclei. 

In Fig. 3, the internal energies, E, are plotted relative 
to the ideal internal energy, Eq. PIMC and DFT-MD 
results for {E — Eq)/Eq agree to better than 0.04 in the 
range of 2.5-7.5x10^ K for water and carbon. Conver- 
gence tests show that results are equally well converged 
in 24-atom and 8-atom simulation cells. PIMC extends 
the equations of state to the weakly interacting plasma 
limit at high temperatures, in agreement with the Debye- 
Hiickcl model. The DFT-MD and PIMC methods to- 
gether form a coherent equation of state over all temper- 
atures. 

Figure 4 shows nuclear pair-correlation functions for 
carbon and water using PIMC and DFT-MD. Fig. 4(a) 
demonstrates the sensitive temperature dependence of 
structural properties for carbon. Water pair-correlations 
are shown in Fig. 4(b) at a single temperature of 7.5 x 10^ 
K. Simulations use a 24-atom simulation cell size, which 
converges pair-correlation curves to better than 10%. 
The PIMC and DFT pair correlation curves essentially 
lie on top of each other with the maximum deviation be- 
ing 17% for carbon at r=0.63 A. The results demonstrate 
that PIMC and DFT predict consistent structural prop- 
erties in addition to the equation of state. 

In conclusion, we have developed an all-electron path 
integral Monte Carlo method using free-particle nodes 
that allows for calculations of materials composed of first- 
row elements and mixtures thereof. Our computations 
of pressures, internal energies, and pair-correlation func- 
tions for water and carbon demonstrate that PIMC and 
DFT can cross- validate each other in a commonly acces- 
sible temperature range and provide an accurate, coher- 
ent equation of state ranging from ambient conditions to 
the plasma limit. The excellent agreement between our 
PIMC method and DFT-MD validates the use of free- 
particle nodes for partially-ionized first-row elements and 
the use of zero-temperature exchange correlation func- 
tional at high temperature. 
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